##############################################################################
# Assembly
##############################################################################
order = 1
phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
D = pde.int.assemble3(MESH, order = order)
# R0, RSS = pde.h1.assembleR3(MESH, space = 'P1', faces = 'ambient_face')
M = phi_H1 @ D @ phi_H1.T
K = dphix_H1 @ D @ dphix_H1.T +\
dphiy_H1 @ D @ dphiy_H1.T +\
dphiz_H1 @ D @ dphiz_H1.T
# Kn = RSS @ K @ RSS.T
phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)
M_Hcurl = phix_Hcurl @ D @ phix_Hcurl.T +\
phiy_Hcurl @ D @ phiy_Hcurl.T +\
phiz_Hcurl @ D @ phiz_Hcurl.T
K_Hcurl = curlphix_Hcurl @ D @ curlphix_Hcurl.T +\
curlphiy_Hcurl @ D @ curlphiy_Hcurl.T +\
curlphiz_Hcurl @ D @ curlphiz_Hcurl.T
C_Hcurl_H1 = phix_Hcurl @ D @ dphix_H1.T +\
phiy_Hcurl @ D @ dphiy_H1.T +\
phiz_Hcurl @ D @ dphiz_H1.T
curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
phix_Hcurl_P0, phiy_Hcurl_P0, phiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = 0)
dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)
KR = R.T@K_Hcurl@R
MR = R.T@M_Hcurl@R
r = jx_hdiv @ D @ phix_Hcurl.T +\
jy_hdiv @ D @ phiy_Hcurl.T +\
jz_hdiv @ D @ phiz_Hcurl.T
# # use this instead if you use the H1-L2 approx for the current J
# r = jx @ D @ phix_Hcurl.T +\
# jy @ D @ phiy_Hcurl.T +\
# jz @ D @ phiz_Hcurl.T
cholKR = chol(KR)
a = cholKR.solve_A(R.T@r)
a = R@a
Bx = curlphix_Hcurl_P0.T @ a
By = curlphiy_Hcurl_P0.T @ a
Bz = curlphiz_Hcurl_P0.T @ a